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ABSTRACT 

We describe an improved, practical method for constructing galaxy models that 
match an arbitrary set of observational constraints, without prior assumptions 
about the phase-space distribution function (DF). Our method is an extension 
of Schwarzschild's orbit superposition technique. As in Schwarzschild's original 
implementation, we compute a representative library of orbits in a given potential. We 
then project each orbit onto the space of observables, consisting of position on the sky 
and line-of-sight velocity, while properly taking into account seeing convolution and 
pixel binning. We find the combination of orbits that produces a dynamical model 
that best fits the observed photometry and kinematics of the galaxy. A new element of 
this work is the ability to predict and match to the data the full line-of-sight velocity 
profile shapes. A dark component (such as a black hole and/or a dark halo) can easily 
be included in the models. 

In an earlier paper (Rix et al.) we described the basic principles, and implemented 
them for the simplest case of spherical geometry. Here we focus on the axisymmetric 
case. We first show how to build galaxy models from individual orbits. This provides a 
method to build models with fully general DFs, without the need for analytic integrals 
of motion. We then discuss a set of alternative building blocks, the two-integral 
and the isotropic components, for which the observable properties can be computed 
analytically. Models built entirely from the two-integral components yield DFs of the 
form f{E,Lz), which depend only on the energy E and angular momentum L^- This 
provides a new method to construct such models. The smoothness of the two-integral 
and isotropic components also makes them convenient to use in conjunction with the 
regular orbits. 

We have tested our method, by using it to reconstruct the properties of a 
two-integral model built with independent software. The test model is reproduced 
satisfactorily, either with the regular orbits, or with the two-integral components. This 
paper mainly deals with the technical aspects of the method, while applications to the 
galaxies M32 and NGC 4342 are described elsewhere (van der Marel et al., Cretton & 
van den Bosch). 



Subject headings: black hole physics — galaxies: elliptical and lenticular, cD 
galaxies: kinematics and dynamics — galaxies: structure. 



1. Introduction 

In order to understand the structure and dynamics of a galaxy, one needs to measure 
tiie total gravitational potential as well as the phase-space distribution function (DF) of the 
constituent stars. The DF specifies the distribution of the stars over position and velocity, and 
hence provides a full description of the galaxy. For a particular galaxy, one needs to explore 
which combinations of potential and DF are consistent with the available observations (surface 
brightness and kinematics). Several methods have been devised to tackle this problem. 

The direct calculation of the DF generally requires analytic knowledge of the integrals of 
motion, and has been restricted in the past to a number of special cases: (i) spherical or other 
integrable potentials (e.g., Dejonghe 1984, 1986; Bishop 1987; Dejonghe & de Zeeuw 1988; Gerhard 
1991; Hunter & de Zeeuw 1992); (ii) nearly integrable systems where perturbation theory can be 
applied (Saaf 1968; Dehnen & Gerhard 1993); or (iii) the subset of axisymmetric models in which 
the DF is assumed to depend only on E and L^ (Hunter & Qian 1993; Dehnen & Gerhard 1994; 
Kuijken 1995; Qian et al. 1995, hereafter Q95; Magorrian 1995; Merritt 1996b, hereafter M96b). 
Numerical calculations of orbits in axisymmetric potentials have shown that most of the orbits 
admit a third integral, which in general is not known analytically (e.g., Ollongren 1962). There 
is no a priori physical reason to expect the DF to depend only on the two classical integrals, and 
in fact, there are indications for both elliptical galaxies (Binney, Davies & Illingworth 1990) and 
halos of spirals (Morrison, Flynn & Freeman 1990) that the DF must depend also on the third 
integral. For the solar neighborhood it has been known for a long time that there must be such a 
dependence (e.g., Binney & Merrifield 1998). 

Schwarzschild (1979, 1982) devised an elegant method to circumvent our ignorance of analytic 
integrals of motion and to build numerically self-consistent equilibrium models of galaxies. 
Richstone (1980, 1984) used this technique to construct axisymmetric scale-free models. It 
was applied to a variety of models (spherical, axisymmetric and triaxial) by Richstone and 
collaborators (see e.g., Richstone &: Tremaine 1984, 1985; Levison & Richstone 1985, 1987; Katz & 
Richstone 1985). Pfenniger (1984) used Schwarzschild's method to build two-dimensional models 
of barred galaxies and Merritt & Fridman (1996) and Merritt (1996a) used it to build a number 
of triaxial models with cusps. Zhao (1996b) modeled the Galactic bar using similar techniques. 
Schwarzschild's original experiment reproduced self-consistently a triaxial mass distribution, but 
as shown by Pfenniger (1984), one can easily include kinematic constraints in the models. Levison 
& Richstone (1985) modeled the observed mean line-of-sight velocities V and velocity dispersions 
a to estimate the amount of counter-rotation in some well-observed galaxies. 

Recent advances in detector technology have made it possible to measure full line-of-sight 
velocity profile (VP) shapes, instead of only the first two moments V and a (e.g., Franx & 
Illingworth 1988; Rix k White 1992; van der Marel k Franx 1993, hereafter vdMF; Kuijken k 
Merrifield 1993). This provides further constraints on the dynamical structure of galaxies. Rix et 
al. (1997, hereafter R97) took advantage of this development, and extended Schwarzschild's scheme 
to model VP shapes. They applied it to spherical models for the EO galaxy NGC 2434, and showed 
that the observations imply the presence of a dark halo. Here we consider axisymmetric models 
and show how to use the extended Schwarzschild method to construct fully general three-integral 



models that can match any set of kinematic constraints. Independent implementations of the 
software were written by N.C. and R.v.d.M. A summary of this development is given by de Zeeuw 
(1997). In an earlier paper (van der Marel et al. 1998, hereafter vdM98; see also van der Marel 
et al. 1997) we applied this modeling technique to the compact E3 elliptical M32, for which 
previous modeling had suggested the presence of a central massive black hole (BH) (e.g., Q95; 
Dehnen 1995). Cretton & van den Bosch (1999) describe an application to the edge-on SO galaxy 
NGC 4342. Other groups are in the process of developing similar techniques to the one described 
here (e.g., Richstone et al. 1997; see also: Emsellem, Dejonghe & Bacon 1999; Matthias & Gerhard 
1999). 

This paper is organized as follows. In Section Q we describe step by step how to construct 
the models (see Figure H). We first discuss the mass models that we consider (Section |2.lD . We 
describe how we choose a grid in integral space that yields a representative library of orbits 
(Section ^^), how these orbits are calculated numerically (Section |2.3D , how their properties are 
stored on a number of grids (Section |2.4| ), and how we model all aspects of the data taking and 
analysis, such as seeing convolution, pixel binning, and extraction of VPs (Section |2.5D . We then 
present the method that we employ to determine the non-negative weight of each orbit (i.e., the 
number of stars traveling on each orbit), such that the global superposition of orbits produces a 
consistent model that best fits the observations (Section \i.6\ ). Lastly, we discuss how we include 
optional smoothness constraints in the models (Section |2.7| ). In Section ^ we describe a set of 
alternative building blocks, the two- integral and isotropic components, for which the observable 
properties can be computed analytically. The smoothness of these components makes them a 
convenient tool to use in conjunction with the regular orbits described in Section ^. Models can 
also be built entirely of these components, to obtain models with DFs of the form f{E,Lz) or 
f{E). In Section ^ we describe the tests that we have performed to establish the accuracy of our 
method. We present our conclusions in Section |5[ 



2. Construction of Dynamical Models 

2.1. Mass Model 

We study dynamical models in which all relevant quantities are axisymmetric, and symmetric 
with respect to the equatorial plane z = 0. It is sufficient to have the total gravitational potential, 
<^ = <1>^ -|- ^I'dark, and the forces, V<1> available and tabulated on a grid, such that their values at 
any point can be recovered through interpolation. This is important, because the structure of real 
galaxies can be very complicated, and is not always well described in terms of analytical functions. 

While the method works for arbitrary radial density profiles, it proves convenient for the 
purpose of presenting and testing our technique to consider models in which the mass density of 
the luminous material, p^, does have an analytical form: 



p,{R,z)=p{s)=po[-\ (l+[-J j (l+[-J j > (1) 

where s is defined as s"^ = R^ + [z/q]'^. This is an axisymmetric generalization of the spherical 



models studied previously by, e.g., Dehnen (1993), Tremaine et al. (1994) and Zhao (1996a), and 
includes the (a, (3) models of Q95 as a limiting case. The model has a constant axial ratio q that 
does not vary with radius. The parameters b and c are characteristic lengths. At small radii 
{r <^ b) the density has a central cusp with logarithmic slope a (when a < 0). At intermediate radii 
{b <^r -^ c) the density falls off as py, oc r'^+iP ^ while at large radii (r ^ c) p^ oc r°'+iP+'^° , When 
viewed at an inclination angle i, the isophotes are ellipses of axial ratio q' = (cos^ i + q^ sin^ i)^'^. 
The luminosity density is j = p*/T, where T is the average mass-to-light-ratio of the luminous 
material, which we assume to be constant. 

For these models, the gravitational potential and the associated radial and vertical forces can 
all be obtained from one-dimensional (usually numerical) integrals (cf. eqs. [2.10]-[2.12] of Q95). 
We calculate the potential and forces in this way and tabulate them on a fine polar (r, 9) grid, 
with logarithmic sampling in radius and linear sampling in the angle. These tabulated values are 
used for the subsequent orbit calculations. It is straightforward to add the contributions from 
a dark component to the potential and the forces, as required for models with, e.g., a BH or a 
dark halo. In the case of a BH these contributions need not be tabulated, because they are known 
analytically. 

For density distributions that are not stratified on similar concentric spheroids one must use 
more general techniques to calculate the gravitational potential and the associated forces. One 
possibility to determine these, while at the same time fitting a complicated surface brightness 
distribution, is to use a Multi-Gaussian Expansion (Emsellem et al. 1994). We do this in our 
modeling of the SO galaxy NGC 4342 (Cretton & van den Bosch 1999). Another possibility is to 
obtain p^ through non-parametric deprojection of an observed surface brightness distribution (e.g., 
Dehnen 1995), and calculate the potential from a multipole or other expansion (e.g., Hernquist & 
Ostriker 1992; Zhao 1996a). 



2.2. Choice of Orbits 

The results obtained with our modeling technique should not depend on the details of the 
orbit library. To achieve this, the library must represent the full variety of orbits in the given 
potential. In this section we describe how we have chosen to select orbits in order to fulfill this 
requirement. 

In axisymmetric models, all orbits conserve at least two isolating integrals of motion: the 
energy E and the vertical component of the angular momentum L^. Numerical studies have shown 
that many orbits conserve an additional third isolating integral Is, which is usually not known 
analytically (see e.g., Ollongren 1962; Innanen & Papp 1977; Richstone 1982; Dehnen & Gerhard 
1993). These regular orbits are specified completely by the integrals of motion, and can be labeled 
by the values of -B, L^ and I^. 

For each energy E^ there is one circular orbit in the equatorial plane, which has radius R^ 
and velocity V^ = Rc{d^/dR)m^fi\. The angular momentum of this orbit, RcVc, is the maximum 
angular momentum at the given energy: L^ax{E). We sample the energies in the model by 



adopting a logarithmic grid in R^. Each Re defines an energy E through the imphcit relation 
E = $(i2c,0) + 2^c- The orbits in the model have Re € [0, oo). However, it is sufficient to adopt 
a grid of A''^; values that covers only a finite range, i?c,min to i?c,maxi chosen so as to contain all 
but a negligible fraction of the total mass of the system. At each energy we sample the range of 
possible Lz values by adopting a grid in the quantity rj = Lz/L^ax{E) {rj e [—1, 1]). Orbits with 
both L^ > and L^ < are included in the library, but the Lz < orbits need not be calculated; 
they are obtained from the L^ > orbits by reversing the velocity vector at each point along the 
orbit. We have calculated orbits for A''^ values of r], spaced linearly between ei and 1 — ei, where ei 
is a small number. For numerical reasons, the special values rj = (radial orbits) and r] = 1 (the 
circular orbit in the equatorial plane) are presumed to be represented by their closest neighbors 
on the grid, but are not included explicitly. 

In an axisymmetric potential the orbit reduces to a two-dimensional motion in the meridional 
{R,z) plane in an effective gravitational potential $eff = ^ + ^L'I/R'^ (e-g-, Binney & Tremaine 
1987, hereafter BT). For fixed {E,Lz), the position of a star is restricted to the region bounded 
by the 'zero- velocity curve' (ZVC), defined as the curve of values {R, z) that satisfy the equation 
E = $eff5 ^-ncl hence vr = Vz = 0. Figure || illustrates ZVCs in the meridional plane. A regular 
orbit admits a third isolating integral, /a, that restricts its motion to a sub-region of the full region 
of phase-space accessible at the given {E,Lz). This is illustrated in Figure ^, which shows a regular 
orbit viewed in the meridional plane. In our method we have chosen a numerical representation of 
Is that can be used to label the orbit. Every orbit with Lz ^ touches the ZVC (OUongren 1962). 
As suggested by Levison & Richstone (1985), we take the R coordinate of the 'turning point' on 
the ZVC (i.e., the intersection of the orbit with the ZVC), denoted by Rzvc, as the third parameter 
to specify the orbit. At every {E,Lz) there is exactly one orbit that touches the ZVC at only one 
value, -Rthini of R: the so-called 'thin tube' orbit (see Figure ^). All other regular orbits touch the 
ZVC for at least two values of Rzvc, one sinaller than -Rthin and one larger than i?thin- To sample 
the orbits at a given (E,Lz), we calculate trajectories that are started with vn = Vz = from 
the ZVC, at a given radius i?zvc (this radius determines v^ according to v^ = Lz/R). Not every 
orbit launched in this way necessarily admits a third integral, since irregular orbits also touch the 
ZVC. Our orbit library therefore includes both regular and irregular orbits, and, as we shall see 



in Section 2.3, we have found it unnecessary to distinguish between them in all our tests and 
applications to date. To reduce redundancy in the library it is sufficient to consider only orbits 
with i?zvc G [Rihm-, Rm&x\i where i?max is the radius at which the ZVC intersects the plane z = d. 

For the orbit library we have chosen to use A''/3 values of -Rzvc- Each point (i^zvc, -^zvc) on 
the ZVC is determined with the help of an angle w, which is sampled linearly between and 
WtYiin (see Figure ^). For numerical reasons, the special values i?zvc = -Rthin (thin tube orbit) and 
^zvc = ^max (equatorial orbit) are presumed to be represented by their respective closest neighbor 
on the grid, but are not included explicitly. Finding the starting point for the periodic orbit, -Rthin; 
is straightforward (see e.g. Pfenniger & Friedli, 1993). 

It is sufficient to calculate only orbits that are started from the ZVC with z > 0. Orbits 
started with z < are obtained from those started with z > by reversing the sign of z and Vz at 
each point along the orbit. Most orbits are themselves symmetric with respect to the equatorial 
plane (see e.g., Figure |3|), so that this operation is redundant. However, this is not true for, e.g.. 



the orbits parented by the 1/1 resonance between the R and z-motion (see Figure 11 below, 
or Figure 8 of Richstone 1982). Since we are only interested in constructing models that are 
symmetric about the equatorial plane, we do not view the orbits started with z > and z < 
from the ZVC as separate building blocks, but instead we consider only their sum. 

The grid in {Rc,r], R^vc) completely specifies the orbit library. Appropriate choices for the 
parameters that characterize this grid are discussed in Section Q. 



2.3. Orbit Calculation 

For each {Rc,r], R^^c) we calculate a trajectory, started from the ZVC as described in 
Section ^^. We have used several standard integration algorithms, including the Bulirsch-Stoer 
integrator (Press et al. 1992) and the Runge-Kutta-Fehlberg algorithm (Fehlberg 1968). We have 
experimented with both and found equivalent results. The former algorithm was used in vdM98. 
Here we use the Runge-Kutta-Fehlberg algorithm. 

The results of the orbit calculations were used to approximate the 'orbital phase-space 
density' for each trajectory. Each phase-point along a calculated orbit was assigned a weight equal 
to the time step at that point, divided by the total integration time. 

This procedure results in density distributions in phase-space, DFtraj and its corresponding 
spatial density ptraj- Orbits were calculated in the meridional plane, but all six phase-space 
coordinates are needed. The azimuthal velocity v^j, = L^/R is completely specified. However, 
for projection onto the sky, the azimuthal angle (p E [0, 2tt] is also required. The distribution of 
stars over (p is homogeneous, because of the assumed symmetry. The weight at each time step 
was therefore divided into a number of equal 'sub- weights', and each was assigned a random 
(j). Furthermore, each sub-weight was divided in two, and one of the two parts was assigned 
phase-coordinates with {z,Vz) multiplied by —1. This corrects for the fact that only orbits started 
with z > from the ZVC were calculated (cf. Section |2.2| ). The trajectories should be integrated 
long enough so that the orbital phase-space densities no longer change significantly with time. 
Pfenniger (1984) proposed to check directly for the convergence of the orbital mass distribution. 
However, this may take a very long time, especially for orbits that are unusually close to a 
high-order resonance, for orbits at very large radii or at very small radii close to a BH, and for 
irregular orbits (see also Merritt &: Fridman 1996). We have used a cruder approach, in which we 
calculated each orbit for a fixed number (~ 200) of characteristic orbital periods. We found this 
to be sufficient for our purpose; longer integrations yield final models that are not significantly 
different. This is because the 'noise' in our modeling is dominated by the representation of phase 
space through a coarse discrete grid. 

Orbits can have sharp edges in both the spatial and velocity dimensions. We found that 
a simple scheme to obtain smoother densities yielded slightly more accurate results for the 
final orbit superposition. To take into account the fact that each energy E in the orbit catalog 
represents all energies in some bin [Ei,E2] around it (defined by the choice of energy grid), 
a random energy E was drawn from the range [Ei,E2] for each normalized time step. The 



corresponding phase-coordinates (r, v) were then translated to the energy E, by replacing them by 
{[Rc/Rc]r, [Vc/Vc]v), where Re, Re and Ve, Ve are the radii and circular velocities of the circular 
orbits at the energies E and E, respectively. This "dithering" approach is only approximately 
correct (it assumes that the potential is locally scale- free), but was found to work well in practice. 



2.4. Storing the Orbital Properties 

For each orbit we store both the intrinsic properties and the projected properties. The 
intrinsic properties are necessary to test for consistency of the final model. We store ptraj on an 
(r, 6) grid in the meridional plane, logarithmic in r and linear in G [0, ^]. Angles G > \ need not 
be stored separately, because of symmetry with respect to the equatorial plane. We also store the 
lowest-order velocity moments of each orbit {ptra.i{vi), ptrai{viVj),i,j = r,9,(p) on the same grid, so 
as to be able to study the intrinsic dynamical structure of the final model. 

The projected properties are necessary for comparison to observable quantities, such as the 
projected surface brightness and line-of-sight VP shapes. Only three coordinates of phase-space 
are available for comparison with observations: the projected positions x' , y' (which we choose to 
be aligned with the photometric major and minor axis), and the line-of-sight velocity, fios(= Vz'). 
Given an inclination angle i of the galaxy (i = 90° means edge-on), these are related to the usual 
cylindrical coordinates (i?, z, 0) in the following way: 

X = RsiiKJ), 
y' = —Rcosicoscj) + zsmi, 
■wios = {vficoscf) — v^s\ii(j))s\iii + VzCosi. (2) 

To have the projected properties of the orbits accessible, we store their phase-space densities 
both on an {r',9') grid on the projected plane of the sky (with similar properties as the intrinsic 
{r,6) grid), and on a Cartesian (x',y',fios) data cube (see Section 2^). The former is used to 



reconstruct the projected surface density of the model. The latter is used to model observed 
kinematical quantities. The spatial grid size (Ax, Ay) of the (x',y',wios) cube is chosen to provide 
2-5 times higher spatial resolution than the pixel size of the available kinematical observations. 
If observations with very different resolution are available for a galaxy (e.g., very high spatial 
resolution HST data in the central arcsec, and lower-resolution ground-based data out to an 
effective radius), it is best to store the data on two or more cubes with different spatial grid sizes 
and extents. During the orbit calculations we then store the phase-space densities simultaneously 
on all {x',y',vios) cubes. Only the x' > half of each cube needs to be stored, because each orbit 
has the same weight at {x' ,y' ,v\os) as at (— x', — y', — fios)- The size Au of the velocity bins on the 
{x\y' ,v\os) cube(s) must be chosen to provide a proper sampling of the observed VPs. In practice 
we use ~ 50-100 bins between (— A'^o-o'maxj-^o-o'max); where (Tmax is the largest observed dispersion, 
and N„ = 4-8. 



2.5. Modeling Observed Kinematical Quantities 

Point-spread-function (PSF) convolution is essential when comparing model quantities with 
observed kinematical quantities in the central regions of galaxies. Seeing convolution correlates 
information in the two spatial dimensions a;',y', but not in fios: 

FconAx'o, y'o, t;ios) = F 0PSF = JJ F{x', y', t;ios) PSF(:e' - x'^, y' - y'o) dx' dy', (3) 

where F is the function to convolve, PSF is the point-spread function, and -Fconv is the result of 
the convolution of F with the PSF. The final model is a linear superposition of the orbits, so the 
{x',y',vios) cubes for each orbit may be individually convolved with the PSF. As in R97, we do 
the convolution for each velocity slice efficiently by multiplications in Fourier space, using Fast 
Fourier Transforms (e.g.. Press et al. 1992). 

Kinematical data is generally obtained either through small, discrete apertures, along a 
number of slit-positions, or may derive from two-dimensional integral field spectroscopy (e.g.. 
Bacon et al. 1995). Any setup of this kind can be simulated by our models, including possible 
spatial binning along a slit. For each observational 'aperture', we choose the {x',y',v\os) cube with 
the most appropriate cell size, convolve it with the relevant PSF, and bin the results spatially over 
the aperture size. This yields a one-dimensional velocity histogram, for each orbit and for each 
observation. Examples of such 'orbital VPs' are shown in Figure ^ 

Kinematical observations provide information on the line-of-sight VPs 

VP(x', y', i;ios) = /// DF dv,, d^ dz', (4) 



at different positions (x',y') on the projected face of a galaxy. In practice, the normalization of 
YF {x' , y' , viqs) is based on the photometric data. It is often useful to parametrize observed VPs by 
a few numbers. A common choice for such a VP parametrization is the Gauss-Hermite expansion. 
We follow the notation of vdMF, in which the VP is represented as 

VP(t;ios) = ^^E^'^'H, (5) 

^ 1=0 

with 

w = {v,,,-V)/a, aiw) = ^e~^'/^. (6) 



/2tt 
The hi are the Gauss-Hermite moments (hereafter GH-moments) defined by 

/oo 
VP(t;ios) a{w) Hi{w) dv\o. (/ = 0, . . . , L). (7) 

-oo 

Each Hi is a Hermite polynomial (see Appendix A of vdMF). The quantities V and a characterize 
the 'weighting function', a{w)Hi, in the integral (^). When describing observations, V and a are 
usually taken to be the velocity and dispersion of the Gaussian that best fits the observed VP. 
With this choice, /ii = /i2 = 0. GH-moments of higher order describe deviations from a Gaussian. 



10 



Only the moments of order L < 6 are generally measured from galaxy spectra, due to the finite 
spectral (and thus velocity) resolution of the observations. 

If we envision galaxies as consisting of orbital building blocks, then the overall VP measured 
for a given observation is just the superposition of the individual orbital VPs. Similarly, the 
observed GH-moments are just a linear superposition of the GH-moments of the individual orbital 
VPs, provided that the observed V and a for the given observation are used in the weighting 
function a{w)Hi{w). Thus, as described in detail in R97, to fit the kinematical observations we 
may restrict ourselves to solving a linear superposition problem for the Gauss-Hermite moments. 
The constraints are then that hi = h2 = 0, and the hi with / > 3 should equal their observed 
values. It must be stressed that this approach is general, and assumes neither that the observed 
VPs are well-described by the lowest-order terms of a GH-series, nor that the orbital VPs are 
well-described by the lowest-order terms of a GH-series. Nonetheless, if a full non-parametric 
estimate of the observed VPs is available there is no need to restrict the analysis to the lowest-order 
GH-moments. Our technique can just as easily fit the individual velocity bins of the observed VPs. 



2.6. Fitting the Constraints 

Constructing a model consists of finding a weighted superposition of the stellar orbits in the 
library that reproduces two sets of constraints: 

• Consistency constraints for the stellar luminosity distribution. The model should 
reproduce the initially assumed luminous stellar density Pi,{R,z) (Section p.l[) , for each 
cell of the meridional {r,0) grid, for each cell of the projected plane {r',6') grid, and for 
each aperture for which there is kinematical data. In theory, it is sufficient to fit only 
the meridional plane masses, because projected densities are then fit automatically. In 
practice this might not exactly be the case, because of discretization. To circumvent this, 
the projected masses may be included as separate constraints. Note that for axisymmetric 
models the projected density does not uniquely specify the intrinsic density (e.g., Rybicki 
1986; Gerhard & Binney 1996). 

• Kinematical constraints. The model should reproduce the observed kinematics of the 
galaxy, including VP shapes. As discussed, we express this as a set of linear constraints on 
the GH-moments of the VPs (R97). 

Finding the orbit superposition that best fits these constraints amounts to solving a linear 
problem, which can be written in matrix notation as B7 = c (R97). The matrix B contains 
the mass that each orbit contributes to each relevant intrinsic or projected grid cell, and the 
GH-moments that each orbit contributes to each kinematical observation. The vector c contains 
the mass predicted by Pi,{R,z) for each relevant intrinsic or projected grid cell, and the observed 
GH-moments for all kinematical observations; the vector 7, which is to be solved for, contains 
the weight of each orbit, i.e., the total mass of stars on each orbit. These weights should be 
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non-negative. Using the terminology introduced in Section 2.5, the basic Schwarzschild equation 
becomes 

-^ 'orbits 

J2 7,VP., = VP,, (8) 

i=i 
with VPjj the individual VP of orbit j at constraint point i, and VPj the observed VP for the 
same constraint point. 

The superposition problem can be expressed as a non-negative least squares (NNLS) fit for 
the above matrix equation. We have used the NNLS routine of Lawson &: Hanson (1974) to solve 
it. The NNLS routine finds a combination of non-negative orbital occupancies (which need not be 
unique) that minimizes the usual C^ norm || B7 — c||. This norm can be viewed as a x^ quantity 
that measures the quality of the fit to the constraints. The NNLS routine always finds a best 
solution. It need not be acceptable in light of the observations; this must be assessed through the 
X^ of the best fit, and by comparison of the model predictions to the constraints. 

As is customary in least-squares fitting, the model predictions for each constraint and the 
actual constraint values (the elements of the vector c) are weighted by the errors in the constraints. 
Observational errors are available for the kinematical constraints. In principle one would like 
the consistency constraints to be fit with machine precision. It turns out that this is generally 
unfeasible, because of discretization. It was found that models with no kinematical constraints 
could at best simultaneously fit both the intrinsic and the projected masses with a fractional 
error of ~ 5 x 10~'^ (when using ~ 1000 orbits). We therefore assigned fractional errors of this 
size to the masses in the consistency constraints. In principle one would like to include also the 
observational surface brightness errors in the analysis. Unfortunately, this requires the exploration 
of a large set of three-dimensional mass densities (that all fit the surface photometry to within the 
errors), which is prohibitively time-consuming. 



2.7. Regularization 

Our orbit superposition models are not generally smooth in integral space, as a result of the 
'ill-conditioned' numerical nature of the NNLS matrix equation that is being solved. There are no 
physical theories that describe exactly how smooth the DF of a stellar system should be, but some 
degree of smoothness should be expected. Our technique can be extended in a straightforward 
manner to yield smooth solutions, by adding linear regularization constraints to the NNLS matrix 
equation (e.g.. Press et al. 1992; Merritt 1993). This has the same effect as the addition of 
'maximum entropy' constraints (Richstone & Tremaine 1988). For linear regularization, each 
regularization constraint must be of the form 

^Sfc,/7fc = 0±A, (9) 

k 

thus providing an extra row to the matrix equation. The 7^ are the orbital weights that make 
up the vector 7, and I is the number of the regularization constraint. The parameter A sets the 
amount of regularization. Models with A ^ cxd have no regularization, while models with A ^ 
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give infinite weight to the regularization constraints. Alternatively, one may view this as adding a 
term A IIS7II to the norm || B7 — c\\ that is minimized by the NNLS routine, where S is the matrix 
with elements {s^^i}, and A = 1/A is a regularization parameter (Zhao 1996b). 

Many choices are possible for the matrix S, with the only requirement that the norm IIS7II 
should provide a measure of the smoothness of the solution. Our choice is based on the fact that 
we consider the NNLS solution 'y{Rc,r], R^^c) to be "smooth" if the second derivatives of the 
(unitless) function ^{Rc,ri,Rzvc)/'yo{Rc) ai'e small. Here the "reference weights" 7o(-Rc) are a rough 
approximation to the energy dependence of the model. These are determined beforehand, e.g., by 
studying the spherical isotropic limit of the given mass density. We view the three-dimensional 
numerical grid in integral space as a Cartesian lattice, and we approximate the second derivatives 
by second order divided differences (eq. 18.5.10 of Press et al. 1992). We assume that the distance 
between adjacent grid points on the lattice is unity, independent of the carthesian direction in 
which they are adjacent. This (arbitrarily) solves the problem that the axes of the grid in integral 
space have different units and yields three regularization constraints for each grid point {i,j,l) 
that is not on a boundary: -7i-ij,i + 27^^,; - Ji+ij^i = 0, -^i,j-i,i + 2-fij^i - 7ij+i,z = 0, and 



3. Two-integral Components and Isotropic Components 

3.1. Definition 

Individual regular orbits correspond to building blocks with a DF proportional to 
6{E - Eq) 5{L-^ - L-^fi) 5{h - h,o), for given (^0,-^2,0, ^3,o)- These are not the only building 
blocks that can be used to construct models. One may also use 'two-integral components', which 
correspond to the DF 

f[Eo,L..o] = ^[i^o,L.,o] SiE - Eo) 5{L, - L,,o), (10) 

or 'isotropic components', which correspond to the DF (cf. Richstone 1982) 



fm^Cm^iE-Eo). (11) 



We choose the normalization coefficients Cieo,l^ 0] ^^'^ ^[Eo] such that the total mass of each 



component is equal to unity; explicit expressions are derived in Appendix A.l. 



The two-integral components are smoother building blocks than the regular orbits, since 
they fill completely the ZVC and do not have the sharp edges of the regular orbits. It is useful 
to view them as a particular combination of all orbits that could be numerically integrated at 
the given {EQ,Lzfl), both regular orbits that fill only a subset of the area enclosed by the ZVC 
(therefore admitting 3 independent integrals of motion) and irregular orbits that occupy a larger 
area (admitting only 2 integrals; note that an irregular orbit does not necessarily fill the entire 
phase-space region defined by (£^01-^^,0); see Merritt & Valluri 1996 for a discussion of the triaxial 
case). Similarly, an isotropic component is a weighted combination of all two-integral components 
(i.e., all orbits) at the given energy Eq. The region in space occupied by such a component is 
bounded by the equipotential surface ^{R, z) = Eq. 
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The two-integral and isotropic components are useful, because their properties can be 
calculated semi-analytically. By using only two-integral components in the NNLS orbit 
superposition, one can construct f{E,Lz) models for arbitrary spheroidal potentials. This provides 
a new and convenient way of constructing such models, which adds to the several techniques 
already in existence for this purpose (Hunter &; Qian 1993; Dehnen & Gerhard 1994; Kuijken 
1995; Magorrian 1995; M96b). Using only isotropic components in the NNLS orbit superposition 
is generally less useful, because these components follow equipotential surfaces, which are rounder 
than isodensity surfaces. Thus, they cannot be used to build self-consistent isotropic axisymmetric 
models. R97 describe how to use them to build spherical isotropic models. Alternatively, the 
two-integral and isotropic components may be used in the superposition in conjunction with the 
regular orbits. This has two advantages. First, these components are smoother, and their inclusion 
therefore reduces numerical noise that arises from the discrete representation of phase space (see 
also Zhao 1996b). Second, addition of these components provides a way to include all irregular 
orbits in the models. 



3.2. Velocity Profiles 

The VP of a two-integral component is obtained by substitution of the DF of equation ([lC|) 
into equation (^). The resulting integral may be written as 

VP[i?o,L.,o] i^'^ y'' ^'los) = C[Eo,L,,o] J d^' J[Eo,L.fi] ' (12) 

(13) 



where 



"" " ' " * ' ■ ' " we 



is the Jacobian for the change of variables from {vx',Vyi) to {E,Lz)- In Appendix A. 2 
give an explicit expression for this Jacobian. The integration in equation ([l2|) extends over 
those z' for which there exist velocities {vx',Vy/) such that E{x' ,y',z',Vx',Vy',vios) = Eq and 
Lz{x' ,y' ,z' ,Vx',Vy',v\os) = Lzfi. Similarly, the VP for an isotropic component may be written as 



VP[s„] (x', y\ 7;ios) = CyE,] J dz' J dL, J^Eo,l.^ ■ (14) 

The projected density for a two-integral or isotropic component, at projected position {x',y'), is 
obtained as the integral of VP(x',y',fios) over vios- 

Equations (12) and (^) can be used to calculate the VPs of the two-integral and isotropic 
components through numerical quadratures, without the need for calculating orbital trajectories. 
The only difficulty lies in finding the domain of integration in z'. We illustrate this for the case of 
an edge-on system {i = 90°). In this case J[Eo^l^ o] = {z'vyi)~^, with 



V = W 2(^0 - $) - -^—} - < (15) 
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(cf. eqs. [ Al3| , [Al^ ). We will refer to the expression under the square-root as W. The integration 



in equation (|l^) extends over those z' for which W > 0. We find the roots of W numerically. We 
start by finding the roots of 2{Eq — <!>) — v^^^ = 0. This gives an interval that encompasses all real 
roots of W, because 2{Eq — $) — vf^^ > W. Then we subdivide this interval in many (~ 100-1000) 
small segments, and check whether the sign of W differs at the ends of each segment. If it does, 
we find the root between these two points through bisection. The continuity of the resulting VP 
was used to check whether all required integration domains in z' were found. For the potentials 
studied here, we typically find two or four roots. 

The VP calculation for edge-on isotropic components is less complicated. The Jacobian is 
quadratic in L^, and can be written as J[Eo,Lz] = K^t ~ ^z){Lz — L~)]~^'^. One can show that L~ 
and Lf are real if 2{Eq — <!>) — vf^^ > 0. In this case the integral over dL^ in equation ([l4 ) extends 
from L^ to Lf, and is always equal to vr. Thus, if there is a ^^^ax ^^^ which 2{Eq — <!>) — vf^^ = 0, 
then 

yP[Eo]ix',y',vios) = 2TrC\Eo]z'^^^. (16) 

If there is no such z^ax (i-^-j if ^los exceeds the escape velocity at the tangent point), then 
YF[Eo]ix',y',vios)=0. 

Figure |5| shows examples of the VP of a two-integral component along the major and minor 
axes of an edge-on system. 



4. Tests 

4.1. The Test Model 

The most useful tests for our axisymmetric implementation are those for which the results can 
be compared to analytical results, or to semi-analytical or numerical results that were obtained 
with an independent algorithm. Models with f{E,Lz) DFs have been widely studied in the 
past five years (e.g., Evans 1993, 1994; Hunter & Qian 1993; Dehnen & Gerhard 1994; Evans 
& de Zeeuw 1994; Kuijken 1995; Q95; Magorrian 1995; M96b). Their properties can be derived 
semi-analytically, and a variety of algorithms and numerical implementations have been presented 
to derive the DF f{E,Lz) that generates a given luminous mass density p-^{R,z) in a given 
potential ^{R,z). These models therefore provide an ideal test case. Here we present two tests 
where we use our method to reproduce the properties of an edge-on f{E,Lz) model. 

We consider a model with a luminous mass density of the form ([l|), with parameters: 
a = -1.435, p = -0.423, 7 = e = 2.0, 6 = -1.298, b = 0.55", c = 102.0", q = 0.73, 
Po = Jo^MQ/L^y, jo = 0.463 x lO^LQy pc~^ (for an assumed distance of 0.7Mpc). We calculate 
the potential of the test model, $ = $^ -)- $dark> under the assumption that T = 2.5, and with the 
option of a nuclear BH ($dark = —GMbh/t) of mass Mbh = 3 x 10^ M©. All these parameters are 
based on the application of our technique to the case of the galaxy M32, which was presented in 
vdM98. This analogy with M32 was chosen mainly to demonstrate the accuracy of our method for 
a realistic galaxy model. 
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The luminous mass density p-^{R,z) and potential ^(R,z) determine uniquely only the 
even part feveniE,Lz) of the DF, f{E,Lz) = feven{E,Lz) + foddiE,Lz). For our test model 
we specify (arbitrarily) the extreme case that foddiE, L^) = fcven{E,Lz) for L^ > 0, and that 
fodd{E,Lz) = -/even(^, ^2) for L^ < (and by definition, /odd(-E',0) = 0). Thus, the f{E,L^) 
test model is 'maximally rotating': all the stars are rotating in the same sense and have L^ > 0. 

First, the unique f{E,Lz) DF of the test model was calculated using the approach described 
in vdM98. We will refer to the resulting DF as DFhq (for Hunter & Qian, 1993). The kinematical 
predictions (VPs, GH-moments, etc.) for the test model DF were subsequently calculated using 
the expressions and software of Q95. The Jeans equations were used as in Cretton & van den 
Bosch (1999) to compute the intrinsic second-order velocity moments (u?) and {v'^) = (vg) in the 



meridional plane. Our tests in Sections [4.2| and ^^ are aimed at assessing how well our algorithm 
can reproduce the test model properties thus calculated with independent methods. This allows 
us to test all key aspects of the orbit model construction, including the sampling of integral space, 
orbit calculation, discreteness effects of the spatial grids, projection into the data cubes, seeing 
convolution, and the NNLS algorithm. Hence, it is no great drawback that our tests are restricted 
to two-intergal models. 



4.2. Reproducing the Test Model with Two-integral Components 

We first describe tests of the extended Schwarzschild technique with only two-integral 



components. We used an {E,Lz) grid as described in Section 2^, with A''^; = 70, -Rc,min = lO""^'^ 
arcsec, -Rc,max = 10^'^ arcsec, A'^^ = 19, and ei = 0.01. Only components with L^ > were included 
in the superposition; the resulting models are therefore by definition maximally rotating with 
a DF of the form f(E,Lz)- The DF is determined uniquely by the mass density. Kinematical 
constraints are therefore not required in the NNLS fit, but only constraints on the consistency of 



the stellar luminosity distribution (see Section 2.6). For these, the polar {r,6) and {r',9') grids 



in the meridional plane and on the projected plane of the sky (see Section |2.4|) were chosen to 
have 16 bins in the radial coordinate between i?c,min and Rc,ma.x, and Ng = Nqi = 7 bins in the 
angular coordinate (a rather modest resolution, but similar tests with finer grids yielded similar 
accuracies). We semi-analytically (eq. [4- 140b] of BT) calculated the isotropic DF f{E) for the 
spherical version of the test model, and used the corresponding masses on our energy grid as 
reference weights for the regularization (see Section |27 



The NNLS algorithm yields the mass on each {E, L^) grid cell, i.e., the integral of dM / dE dL^ 
over the grid cell. It does not directly yield the DF f{E,Lz), which by definition is the density 
in the six dimensional phase-space. However, for a two-integral model there is a simple relation 
between dM / dEdL^ and the DF f{E,Lz), as derived in Appendix p.l|. With equation (B3) the 



NNLS fit provides an estimate of the DF, which we will denote DFnnls- We compare it with 
DFhq on the same grid, but to avoid possible border effects, we restrict the comparison to the 
Ne = 50 energy grid points with Re between 10~^ and 10^ arcsec (see Figures ^ and for the test 
models with and without a central BH, respectively). 



16 



The DFs agree well over 10 and 20 orders of magnitude, respectively. The inserts show 
the percentage errors in the DF calculation. Note that the largest errors occur at grid points 
that carry little mass, e.g., at large radii. The orbit library in these figures is numbered as 
follows. For each value of Rc{E)j of the energy, L^ runs monotonically from Lz^ram to L^^max, 
covering Mr] (= 19 here) components. The next orbit corresponds to Lz,min of Rc{E)j j^i^ etc. 
This choice of numbering causes the jagged appearance of the DF. Q95 plotted f{E,Lz = 0) 
and f(E,Lz = -^2, max) as a function of E (see their Figure 8), and such curves would appear as 
envelopes in Figures and 1^. 

To assess the influence of smoothing on the accuracy with which our technique recovers the 
DF, we have studied the dependence of the RMS logarithmic residual 



RMS 



logDF 



1 



Ne ^v 



NeN, fr[ f^. 



5]5^(logDFNNLS-logDFHQ)' , (17) 



1/2 



on the regularization parameter A (see Section p.?]) . When A tends to zero, the regularization 
constraints receive infinite weight. This yields a very smooth DF, but one that doesn't fit the 
consistency (mass) constraints very well, and therefore doesn't approximate DFhq very well. At 
the other extreme, when A is very large there is hardly any smoothing, and the mathematically 
ill-conditioned nature of the problem yields a very jagged solution that also doesn't match DFhq. 
Figure ^ shows RMSiogDF(^) for the following three cases: (a) the case in which only the masses 
on the meridional plane (r, 9) grid are included as constraints in the NNLS fit; (b) the case in 
which only the masses on the projected plane (r', 6') grid are included as constraints; and (c) the 
case in which both are included as constraints. In principle, deprojection of the projected mass 
density is unique for an edge-on system, so these approaches should recover equivalent results. 
However, this is not true in practice because of discretization effects. In all three cases, RMSiogOF 
has a minimum near logA ~ 2.0, which is thus the optimum smoothing. The value of RMSiogDF 
at the minimum is only mildly different for the different cases, but (a) yields the slightly better 
results. We have therefore adopted case (a) for all our further test calculations. The minimum 
RMSiogDF is ~ 0.02; this corresponds to a 5% RMS difference between DFnnls and DFhq. The 
mass-weigthed RMS difference. 



RMSdf 



DF 



HQ 



DFnnLS-DFhQ12 3^ 3^ / ffy.^ j3-j3- 
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(18) 



for the model with the optimum smoothing is also ~ 5%. This level of accuracy in the 
determination of the DF is similar to that obtained with other techniques (see Gerhard et al. 1998; 
Matthias & Gerhard 1999). 

Figure ^ compares the predictions for the meridional plane velocity moments to the results 
of the Jeans equations, for the model without a BH and with the optimum smoothing. The top 
and bottom row of the figure show (v'i^)^''^ and (t^^)^ = (^0)^ 1 respectively. We plot separately 
each angular sector of the polar grid in the meridional plane. In each row, the first panel is 
closest to the symmetry axis and the last one is closest to the equatorial plane. Full lines show 
predictions of the extended Schwarzschild technique, and dashed lines the results obtained from 
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the Jeans equations. The model predictions were interpolated between the {E,Lz) grid points to 
get smoother results. Overall the agreement is very good, and better than 1%. This is better than 
the ~ 5% agreement in the DF, because the velocity moments are integrals over the DF (such that 
errors tend to cancel). The errors in the velocity moments are largest near the symmetry axis, 
since in the extended Schwarzschild technique only a few components with very low Lz can reach 
this region of the meridional plane. However, the errors are always < 2kms^^. 



4.3. Reproducing the Test Model with Regular orbits 

The next step in our testing procedure is to try to reproduce the properties of the test model 
with regular orbits, rather than two-integral components. The first obvious question is whether we 
can give the orbit superposition algorithm constraints that force it to generate a model with a DF 
of the form f{E,Lz), which can then be compared to the distribution function DFhq calculated as 
in Section [4.1| . Unfortunately, there is no set of simple linear kinematic constraints that force the 
NNLS algorithm to produce an f{E, L^) model. One can certainly impose the necessary conditions 
that (v^) = {vq) and {vrVg) = 0, but these conditions are not sufficient; an f{E,Lz) model is fully 
determined only by constraints on all its higher order velocity moments (e.g., Magorrian &; Binney 
1994). 

We therefore restrict ourselves here to a simpler test. We calculate an orbit library in the 
gravitational potential of the test model, but do not do a subsequent NNLS fit. Instead, we 
fix the orbital weights ■~fj to those appropriate for an f{E,Lz) model, and merely calculate the 
projected kinematical quantities for some observational test setup, given these orbital weights. 
The results are compared to the same quantities but now calculated from DFhq as described in 



Section 4.1. This tests all of the important parts of our method that were not already tested by 



the calculations in Section 4.2, namely the orbit calculation, the projection into data cubes and 



VPs, and the seeing convolution. 

The main difficulty with this test is that it requires knowledge of the orbital weights for an 
f{E,Lz) model, i.e., of the differential mass density dM / dEdLzdRzvc on the grid of quantities 



{E,Lz,Rzvc) that we use to sample orbit space (cf. Section ZJ). This is not as straightforward as 



the calculation of dM / dE dLz described in Appendix B.l. In fact, the orbital weights can only 
be easily calculated if an explicit expression exists for the third integral, which is not the case for 
our test model. However, if the model has a central BH then at small radii, or high energies, the 
potential is Keplerian and spherical (as it is at large radii, or low energies, because of the finite 
mass of the model). In this potential, all the integrals of motion are known, and these limits are 
therefore analytically tractable. 

At high energies (small radii) the test model reduces to a scale-free axisymmetric mass density 
cusp with an f{E,Lz) DF in a spherical Kepler potential. This limit was studied analytically by 
de Bruijne, van der Marel & de Zeeuw (1996). In this limit, the normalized distribution of mass 
over {r],Rzvc/Rc) at fixed energy, which we will denote as h{r] , R^^c / Re) , is a known function that 
is independent of energy. An explicit expression for h{rj,Rzvc/Rc) is derived in Appendix |B.2| . We 
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use this result to approximate the differential mass density of our test model, restricting ourselves 
to the case with a 3 x 10^ M© BH. First, we calculate the differential mass density of the model 
over energy alone: G{E) = dM / dE = J[dM / dE dL^] dL^; where dM / dEdL^ is obtained from 
equation (p^). Then, we assume that the distribution of mass over {Lz,Rzvc) at fixed energy is 
always the same as in the high-energy limit, so that 

'^^ --G{E)h{ri,R,,jR,). (19) 



dEd7]d{Rz,^/R^ 



This relation is only correct at asymptotically high energies. We found it was sufficiently good at 
energies with Rc{E) ^ 0.5", which is where the mass density of the model is a pure power law and 
where the potential is Keplerian. 

For our test we picked an observational test-setup with the same set of 8 square apertures, 
roughly aligned on the major axis, as used by van der Marel et al. (1997b) in their HST observations 
of M32 (see their Figure 3). These apertures all lie at projected radii R < 0.5", so most of the 
light seen in these apertures originates from stars with energies for which our approximation of 
the differential mass distribution is adequate. We chose the same PSF as in van der Marel et 
al. (1997b), which is a sum of three Gaussians that approximates the HST PSF. Subsequently, 
we picked a grid in {E, Lz, Rzvc) space with Ne = 20, Nyj = 7, and Nj^ = 7, and we calculated 
an orbit library for this grid. Then finally we calculated orbital weights from equation (^), by 
integrating at each point of our {E,Lz, Rzvc) grid the approximation dM/ [dE dr] d{Rzvc/ Re)] over 
the corresponding grid cell. Predictions for the projected kinematics then follow by superposing 
the VPs for individual orbits in the library as in equation (||). 

Figure ^ shows the results thus obtained for the kinematical quantities V, o", /i3 , . . . , /ig • As 
mentioned in Section |^, independent software implementations of the extended Schwarzschild 
method were programmed by both N.C. and R.v.d.M. Dotted curves in the figure show the 
results from N.C.'s software, while dashed curves show the results from R.v.d.M. 's software. 
For comparison, solid curves show the results obtained by direct integration over the known 



DFhq, using the (completely independent) software of Q95 as described in Section iA. The 
RMS difference between the different predictions is ~ 2kms~^ in V and a, and ~ 0.01 in the 
Gauss-Hermite moments. Kinematical data typically have larger observational errors than this, so 
the numerical accuracy of our method is entirely adequate for modeling real galaxies. 



Finally, let us say a few words about the orbits in the library for this test model. Figure 11 



shows the orbits in the library at the energy corresponding to Re = 0.25". Figure 12 shows the 
orbits at the same Re, in the same model, but now without a BH. The orbits are all regular and 
have a stable periodic parent. The parents were determined using surfaces-of-section (see also: 
Richstone 1982; Lees &: Schwarzschild 1992; Evans 1994; Evans, Hafner & de Zeeuw 1997) and 



are indicated in the figures. Most orbits in Figures 11 and |1^ are tubes, and are parented by the 
thin tube. Orbits that are not tubes are indicated. A minority of the low \rj\ orbits in Figure 12 
is parented by higher-order resonances, such as the 3:2 and 4:3 (being the ratio of the R- and 
z-frequencies of the parent). By contrast, most of the low \r]\ orbits in Figure ^ is parented by 
the 1:1 resonance. Thus the orbital structure of the models with and without BH is very clearly 
different. An analysis of the orbital structure of these models as a function of the BH mass is 
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beyond the scope of the present paper, but does seem worth further study. 



5. Concluding Remarks 

In this paper we have described an extension of Schwarzschild's method for building 
anisotropic axisymmetric dynamical models of galaxies. We compute a set of orbits in a given 
mass model and find the non-negative superposition of these orbits that best reproduces a set of 
(photometric and kinematic) constraints. Our method includes the full VP shape as kinematic 
constraint. We parametrize the VP using a GH expansion so that it is specified by a few numbers. 
The modeling method is valid for any kind of parametric (or non-parametric) VP representation 
and properly takes into account the observational setup (seeing convolution, pixel binning, error on 
each constraint). We obtain smooth models by imposing a regularization scheme in integral-space. 

R97 have described in detail several aspects of this method and applied it to the spherical 
case. However, it is not restricted to this simple geometry, and we have described here the 
axisymmetric extension. The mass model used to compute the orbit library may be complex: it 
can have a central density cusp, a stellar disk, a central black hole or an extended dark halo. 
Applications of our code to the flattened systems M32 and NGC 4342 were presented in vdM98 
and Cretton & van den Bosch (1999). 

We have also devised a new semi-analytic method for constructing simpler dynamical models, 
for which the DF has the special form DF = f{E,Lz) or DF = f{E). These DFs are obtained 
by using NNLS with analytic building blocks for which the VPs are obtained by one-dimensional 
quadratures. This technique is general and does not require the density to be expressed analytically 
as a function of the potential, but can be used with any complex mass model. Previous techniques 
assuming DF = f{E,Lz) that are also free of this condition include those of Hunter &: Qian 
(1993), Dehnen & Gerhard (1994), Kuijken (1995), Magorrian (1995) and M96b. 

We have tested our new method by having it reproduce the properties of f{E,Lz) models 
for which the DF and projected properties can be calculated with independent algorithms. This 
allowed us to test all aspects of the superposition method, including the sampling of integral space, 
orbit calculation, discreteness effects of the spatial grids, projection into the data cubes, seeing 
convolution, and the NNLS algorithm. Tests with only two-intergal components reproduced the 
DF with a mass- weighted RMS accuracy of < 5%, and the meridional plane velocity moments 
to better than 2kms^"'^. Tests with a regular orbit library indicated accuracies in the projected 
quantities of ~ 2kms^^ in V and a, and 0.01 in the GH-moments. All the tests that we have done 
indicate that the accuracy of our method is adequate for the interpretation of kinematical data 
obtained with realistic setups. 

Our technique can be extended to triaxiality. Several parts of the method will be unaltered 
for this geometry: the use of projected quantities, e.g., YF (x' , y' , vio^) , the fitting procedure, the 
seeing convolution, etc. However, the orbital structure is much richer than in the axisymmetric 
case. New orbit families appear (e.g., box orbits) as well as numerous chaotic regions associated 
with resonances. During the numerical integration of a trajectory in such a mass model, all six 
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phase-space coordinates need to be computed, since there is no azimuthal symmetry. Consequently, 
the computing overhead is significantly higher for triaxial geometries. Work along these lines is in 
progress. 
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A. Construction of f{E,Lz)- and /(£')-coniponents 
A.l. Normalization 

To determine the normalization coefficients Cieq,Lzo] ^^^ ^[Eo] ™ ^^^ definitions of the two- 



integral and isotropic components (eqs. [|l^, 11|), we seek expressions for the total mass of a single 
component. The phase-space volume in cylindrical coordinates is d'^x d^w = R^ dRd(j)dz dRd(j)dz. 
We use 

R^ = 2[E-'^{R,z)]-^-i'', ^ = LJR\ (Al) 

to switch, at fixed (-R, z), from the variables (i?, i) to {E, L^)- We then find for the total mass of a 
two-integral component: 

rr^[Eo,L.,]^Jd'xd'vf^^^^,^^^=jR'dRjdc^JdzJ^J^Jdzf^^^^^^^^^^ (A2) 

Upon substitution of /r^ ^ ^ from equation (|lO|), the integration over cj), E and L^ becomes 
trivial, and we obtain 

^lEo,L^,o] = ^^C[E„L^^^] JJdRdzJ{2 [Eo - HR,z)] - ^ - i'}"'^' di. (A3) 

The integral over di extends over the region \z\ < imax, where imax is defined as the root of the 
expression in the square root. Therefore, 

^lEo,L.,o] = ^^C^EoMfi] dRdz = 47r q£o,L,,o] // d-R dz, (A4) 

where the remaining double integral is over the region for which Eq — ^{R, z) — (L^ q/2R^) > 0. 
This is exactly the region ^eff.ol-^i-^) < Eq bounded by the ZVC at the given {Eq^Lz^), where 
$eff,o is the effective gravitational potential at the given Lzfl. To obtain imyEo^L^o] = 1) we choose 



Q 



[Eo,L,fi] 



-1 



27r^ <h [Rdz-zdR) , (A5) 



-1 



4^2 dRdz 

'S>cB.oiR'^)<Eo ZVC[£o,iz,o] 

where the second equality was obtained with a variant of Stokes' theorem. 

Following similar arguments, we obtain for the mass of an isotropic component: 

miEo] = Jd'xd'vf^Eo] (A6) 

= STrq^;^] J J dR dz J dLz j'j^^ /^_ .^ = 47r2q^„] jjj dR dz dL„ 

where the effective gravitational potential $cff(^) z) is a function of L^, at fixed {R,z). Evaluation 
of the integral over dL^ yields 

m^Eo] = 8V2 7r^ C[Eo] j R^Eo-^{R,z) dRdz. (A7) 

<i>{R,z)<Eo 
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To obtain m^Eo] = 1, we choose 

C[Eo] = [8V2 t:^ I R^Eq-^{R,z) dRdzy^. (A8) 

'i>{R,z)<Eo 
For the special case of a spherical potential, ^ = '&(r), we have 

m[s„] = 16^/2 ^2q^^^ j r^^E^-^r) dr, (A9) 

<I>(r)<£o 

which can be recognized as the 'density-of-states' function for an isotropic spherical system 
(BT). Calculations for spherical components with DFs proportional to 5{E — Eq) 6{L — Lq) were 
presented in Appendix A of R97. 



A. 2. Velocity profiles 

We derive here the Jacobian J for the transformation from {vxi,Vyf) to {E,Lz), which enters 
into the expressions for the VPs of two-integral and isotropic components (eqs. [|12| , |l^ ]). The 
energy is 

E = i(t;2, + ^2, + yfj + $(^'^ y'^ /). (AlO) 

For inclination angle i, 

X = —y' cosi + z' sini, z = y' sini + z' cosi, (All) 

and the angular momentum is therefore 

Lz = xvy — yvx = Vx'{—y' cosi + z'sini) + Vy'{x' cosi) — t'ios(x'sini). (A12) 

This yields for the Jacobian 



1 
J ' 



x'vx' cos i + y'vy' cos i — z'vyi sin i 



in which Vx' and Vy' are functions of E and Lz determined by: 
Lz — Vy/x' cos i + f los^;' sin i 



(A13) 



Vx' = - '- . /.\ , vp = 2{E-<^)-vi,-vL. (A14) 

—y' cos I + z' sm i ^ 



Substitution of Vx' in the expression for v^, yields a quadratic equation for Vyt : 

'J^' + bVy 



avl, + bvyi + c = 0, (A15) 



where 

a = {—y' cosi + z sm.i) +{x'cosi) , 

h = — 2{Lz +v\osX sin i) X cosi, 

c = VY^^{—y' cosi + z' sini)'^ + {Lz + v\oiiX' suiiY — 2{E — (^){—y' cosi + z' suii)"^ . (A16) 

Therefore, 

2{Lz + v\osx' sin i) x' cos i ± vA 
2[{—y' cos i + z' sin i)^ + (x' cos i)^] ' 
with A = b'^ - 4ac. Equations ( |A13D , ( |A14| ) and ( |A17| ) define the Jacobian J. 



V = or/ ../^„.. , ./„.^.N2 , /^/„„„.^21 ' (A17) 
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B. Relation Between Orbital Weights and the DF 

B.l. dM /dEdL-, for an f{E,L^) model 

For a two-integral model there is a simple relation between the differential mass density 
dM / dEdLz and the DF f{E,Lz). To derive this relation (see also Vandervoort 1984) we write 
the trivial identity 

f{E, L,) = JJ fiEo, L,,o) {f[Eo,L.,o]/ClEo,L.,o]) d^o dL,,o, (Bl) 

where /r^ ^ i and C^Eq^l^ q] are as defined in equation (|^). The total mass of the system is 

M ^ J d'xd'v f{E,L,) = JJ f{Eo,L,^o) im[E,,L.jqEo,L^o]) dEo dL,,o, (B2) 

where the second identity follows upon substitution of equation ( |B1| ) , exchange of the integration 
order, and use of the definition of fn^Eo^L^o] from equation (A2). Substitution of equations (A4) 



and (A5), relabeling of the integration variables from {E(j,Lzfi) to {E,Lz), and differentiation 
then yields 

dM /dEdL, = f{E,L,) X [27r2 / {Rdz-zdR)\. (B3) 

zyc[E,L^] 
The mass weight 7^ for a cell in integral space is 



B.2. Scale-free Density in a Kepler Potential 

We summarize here the asymptotic case of a scale-free axisymmetric density in a spherical 
Kepler potential, which was discussed in detail by de Bruijne, van der Marel & de Zeeuw (1996). 
We adopt the same units as in that paper. In those units, /O^ = s~^ and ^ = 1/r, where s is 
defined as in eq. (0) by s^ = B? + {z/q)^ , with q the axial ratio. The associated eccentricity is 
e = \/l — q^. It is convenient to work with the integrals of motion 

£ = -E, 7?2^l2/lLx(^), C'^^V^Lx(^), (B5) 

where £ is the binding energy, rf G [0, 1] and C^ G [rf, 1]. The unique (even) two- integral DF is 

DF(£:,r/2) = Co5(f)i(eV), (B6) 

where 



2vr^(e-i,|) 



Co = ^_J , ,, , 9{£)=£'~l j(eV) = 3^2(1, ^,^;i,^;eV). (B7) 
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The special functions B and 3F2 are the beta-function and a generahzed hyper geometric function, 
respectively. Upon substitution of ^ = —a this yields the Re ^ limit of our test model 



(Section 4.1); upon substitution of ^ = —a — (P'j) — (Se), it yields the Re ^ 00 limit. 



We wish to calculate the mass weight 7j contained in a cell number j of integral space (see 
eq. P). According to equations (35) and (37) of de Bruijne et al., this is given by: 

^,=llld£dv'dewi£,rj\e)^F{£,rj% w{£,rj^e) = ^ £-'/' iv^y'^' iC'r'^' . (B8) 

cellj 

In the Kepler potential the binding energy is related to the circular radius according to 
£ = l/{2Rc). The ZVC at a given [RcTff) is therefore defined by 

1 = ^ - (^f- m 

A particle on the ZVC has Vr = vq = 0, Lp' = r'^v'^ and L^ = R^vi^. Therefore, we have 
(^2 — (^rrj/R)"^. Combined with the expression for the ZVC this yields 

2_r 2n{R,,jRe) i2 



which is a one-to-one relation if (Rzvc/Rc) is chosen between r] and 1 + ^/T—rf. Substitution 
of the DF from equation (P6D into equation (P8|) , and transformation to the variables logi?c; 



rj e [0, 1] and {R,,c/Rc) G [r/, 1 + y/l^^] yields 

d(iogi^.)d,d(i^...w =^^°- ^t) ''"(^^Hl^^^TTiWiOT/- ^^''^ 

Hence, the normalized distribution of mass over {r],Rzvc/Rc) at fixed energy, which we will denote 
as h{rj,R2vc/Rc), is independent of energy. In particular: 
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EXTENDED SCHWARZSCHILD METHOD 
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Fig. 1. — Flowchart of the extended Schwarzschild method. We find the non-negative superposition 
of the orbits with a least squares algorithm (NNLS). This combination of orbits reproduces a set 
of photometric (surface brightness distribution) and kinematic constraints (VPs). 
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N 




Fig. 2. — Zero Velocity Curves are plotted for 7 values of L^ uniformly sampled between 0.05L^^^ 
and 0.95L^^^ at an energy corresponding to a circular orbit radius Re = 2.0 (indicated by the dot). 



The mass model is the one of our test model of Section 4.1, with no BH. 
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Fig. 3. — A regular orbit, the thin-tube periodic orbit, and the ZVC around them in the meridional 
{R, z) plane for the same mass model as in Figure §. The radii Rmin, Rmax and Rthin are indicated 
(see text). The circular orbit is represented by a black dot. The different starting points on the 
ZVC are shown with open dots and the corresponding angles w and t^thin are indicated. 
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Fig. 4. — VPs as a function of line-of-sight velocity wios (in kms~^) for the two orbits of Figure]^. 
The regular and thin orbits are shown in the top and bottom panels, respectively, for viewing 
through (l"-square) cells along the major axis (left) and minor axis (right), respectively. The 
orbits were not convolved with a PSF. For each panel, different lines correspond to cells at different 
distances r from the center: the full line corresponds to the central cell (r = 0"), the dotted line to 
r = l" , the short dashed line to r = 2" and the dot-short dashed line to r = 3". 
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Fig. 5. — VPs of an individual two-integral component with the same {E,Lz) as the orbits in 
Figure S viewed along the major axis (left) and minor axis (right), respectively. Line types are 
the same as in Figure Q. However, unlike in Figure |^, these VPs were calculated for a point along 
either axis, and were not 'integrated' over cells. 



33 



o 

to 



o 



- 



-10 - 




15 



200 400 600 800 
Component number 



500 

Component number 



\ 

1000 



Fig. 6.— The DF f{E, L^) for the test model with a 3 x 10^ M© BH described in Section |4|. The 
DF from the semi-analytical HQ algorithm and from the extended Schwarzschild technique, using 
only two-integral components, are plotted as function of the component number. The components 
run in order or energy, and in order of L^ for each energy; this causes the jagged appearance of 
the curves. The two curves mostly overlap in the comparison interval (dashed lines). The insert 
shows the relative difference (in per cent) between the two DFs. The agreement between the two 
methods of calculating the DF is acceptable. 
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Fig. 7. — Similar to Figure y, but for the same model without a central BH. 
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Fig. 8. — The RMS logarithmic residual RMSiogOF for the test model of Figure 0, as function of 
the logarithm of the regularization parameter A. The residual measures the difference between 
the DF as calculated with the semi-analytical HQ algorithm, and as calculated with the extended 
Schwarzschild technique, using only two-integral components. The different line-types indicate the 
cases in which only masses in the meridional plane are included as consistency constraints (dotted 
line), in which only projected plane masses are included (dashed line), or in which both are included 
(full line with dots). 
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Fig. 9. — Comparison of meridional plane velocity moments, calculated either with the extended 
Schwarzschild technique using only two- integral components (full curves), or with the Jeans 
equations (dashed curves), for the model of Figure M. The top and bottom row of panels show 
(^i)"^ and {vl)^'"^ = {vg)^'"^, respectively. The meridional {r,6) plane is divided in 7 sectors. In 
each row, the first panel is the sector closest to the symmetry axis and the last panel is the sector 
closest to the equatorial plane. The discrepancies are largest near the symmetry axis, but are 
acceptable everywhere. 



37 




150 



100 



b 



2 4 6 8 

aperture number 




2 4 6 

aperture number 



0.05 



-0.05 




0.1 



0.05 



^ 



-0.05 



-0.1 



0.1 



0.05 




2 4 6 8 

aperture number 



-0.05 



-0.1 




2 4 

aperture number 



2 4 6 

aperture number 



2 4 6 8 

aperture number 



Fig. 10. — Kinematical predictions for the edge-on f{E,Lz) test model with a 3 x 10^ Mq BH 



discussed in Section 43. The kinematical apertures are the same as for the HST observations 
of M32 by van der Marel et al (1997b). They are aligned along the major axis. Data points 
are plotted equidistantly along the abscissa. Dotted and dashed curves are predictions obtained 
with the extended Schwarzschild technique, using the software of Cretton and van der Marel, 
respectively. The solid curves show predictions obtained from direct integration over the DF, using 
the software of Q95. The three curves agree well, demonstrating the numerical accuracy of the 
orbit superposition technique. 



77 



V 






1 \''--"\ 




m 

r---"i 1 




1 r'---"i 1 




1 .-^ 

i; , ■ ~ 
1'— -I 




1 i.-wi 

ll 

1 ]'■■-■'] 




1 1 ..-J 1 

i 
1 I'-— -1 1 




1 1,--. 1 1 

1 l'--"i 1 




-10 12 




-10 12 




-10 12 




-10 12 




-1 D 1 2 




-10 12 




-10 12 





III 
1 1 ""'1 




III 

1 ""1 1 




1 1 1 1 

1 ""'T 




II 

1 ■. -■ ~ 
1 """'1 




III 

\ 
1 1 """1 




7W\ 




1 1 1 1 

-Is 

1 1 """1 1 




-10 12 




-10 12 




-1 1 E 




-10 12 




-10 12 




-1 1 E 




-10 12 





III 
III 




III 

iT^/" 
III 




II 
II 




II 

1 '. 

1 •-_.-■■ ~ 
II 




III 

Hi: 

III 




1 1 1 1 

Hi: 

1 1 1 1 




1 1 1 1 

m 

1 1 1 1 




-1 1 H 




-1 1 H 




-10 13 




-1 1 a 




-10 13 




-10 13 




-1 1 s 





III 

2:1 / J 

ill 




III 

3:1 / Ji 

III 




II 

2:1 /' ^ 




II 

II 




III 

Hi 

III 




1 1 1 1 

Ht: 

1 1 1 1 




1 1 1 1 

Hf: 




-10 12 




-10 12 




-10 13 




-10 12 




-10 13 




-10 13 




-10 13 





III 

2:1 ,■'■""""■. 

Ill 




III 

3:1 1 /""^ 
III 




II 

2:1 /"^ 

II 




II 
1 

— i"^B"" 
1 '-.^ 
1 
II 




III 

1 
III 




1 1 1 1 

%; 

1 
1 1 1 1 




1 1 1 1 

H# 

1 1 1 1 




-10 12 




-10 12 




-10 13 




-10 12 




-10 13 




-10 13 




-10 12 





III 

2:1 /--., 

Ill 




III 

2:1 /■ -.. 
Ill 




II 

II 




II 
1 

1 ■-■■ 
1 

II 




III 

III 




1 1 1 1 

1 
1 1 1 1 




1 1 1 1 

H^: 

1 1 1 1 




-10 12 




-1 D 1 2 




-1 D 1 2 




-10 12 




-1 D 1 2 




-10 12 




-10 12 





III 

2:1 ,-.. 

Ill 




III 

2:1 ,-, 

III 




II 

II 




II 

1 
1 

II 




III 

"2:1 I ^ " 

1 
III 




1 1 1 1 

1 

1 

1 1 1 1 




1 1 1 1 

1 1 1 1 



-10 1 S -1 1 



-1 1 



-1 1 E -1 1 



-> 



w 



Fig. 11. — Examples of orbits at the energy with Re = 0.25", in the test model of Section 4.1 for 
the case with a 3 x 10^ M© BH. The axes for each orbit are in units of Re- Each line corresponds 
to a different value of \r]\ and each column to a different value of the third integral. The ratio of 
the R- and ^-frequencies for the parent orbit is indicated for those orbits that are not parented by 
the thin tube. 
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Fig. 12. — Similar as Figure 11, but now for the same model without a central BH. 



